#include <fstream>
#include "spline.h"
#include <iostream>
#include "json.h"
#include "jsoncpp.cpp"
using namespace std;

int main()
{


  ofstream proC;
  proC.open("proC.txt", std::ios::out);
  
  vector<double> x1;
  vector<double> y1;
  int n1=11;
  for (int i=1;i<=11;i++)
    {
      double x;
      x=-6+i;
      double y;
      y=1/(1+x*x);
      x1.push_back(i);
      y1.push_back(y);
    }

  Bspline fB1(11,x1,y1);
  vector<double> plot1;
  double sub=0.1;
  double f_5=10/pow(26,2);
  vector<double> xd;

  Json::Reader reader;
  Json::Value value;
  
  ifstream fin("CD.json",std::ios::binary);
  if (!fin.is_open())
    {
      cout << "Can not open CD.json!" << endl;
    }
  if(reader.parse(fin,value))
    {
      for (int i=0;i<7;i++)
	{
	  double xt = value["xd"][i].asDouble();
	  xd.push_back(xt);
	}
    }
  fin.close();


  vector<double> SD1(7);
  vector<double> SD2(7);
  vector<double> SD3(7);
  vector<double> SD4(7);
  vector<double> SD5(7);
  vector<double> yd(7);
  vector<double> Es1;
  vector<double> Es2;
  vector<double> Es3;
  vector<double> Es4;
  vector<double> Es5;
  vector<double> error1;
  vector<double> error2,error3,error4,error5;
  
  //cout << "yd---" << endl;
  for (int i=0;i<7;i++)
    {
      yd[i]=1/(1+xd[i]*xd[i]);
      // cout << yd[i] <<" "; 
    }
  //cout << endl;
 
  fB1.answer=fB1.S_23(1,f_5,-f_5);
  proC << "-----" << "cubic Bspline D1" << "-----" << endl;
  
  for (int k=0;k<10;k++)
    {
      double coe_0=216*fB1.answer[k].coe[3]+36*fB1.answer[k].coe[2]+6*fB1.answer[k].coe[1]+fB1.answer[k].coe[0];
      double coe_1=108*fB1.answer[k].coe[3]+12*fB1.answer[k].coe[2]+fB1.answer[k].coe[1];
      double coe_2=18*fB1.answer[k].coe[3]+fB1.answer[k].coe[2];
      double coe_3=fB1.answer[k].coe[3];

      /*
      double x=fB1.answer[k].begin-6;
      double y=1/(1+x*x);
      cout << y <<" ";
      double yy=coe_0+coe_1*x+coe_2*x*x+coe_3*pow(x,3);
      cout << "** " << yy << " "<<endl;
      */

      for(int i=0;i<7;i++)
	{
	  if (fB1.answer[k].begin-6<xd[i] && xd[i]<=fB1.answer[k].end-6)
	    {
	      SD1[i]=coe_0+coe_1*xd[i]+coe_2*xd[i]*xd[i]+coe_3*pow(xd[i],3);
	    }
	}
  
      for(int s=0;s<10;s++)
	{
	  double plotx=fB1.answer[k].begin-6+s*sub;
	  double ploty=coe_0+coe_1*plotx+coe_2*plotx*plotx+coe_3*pow(plotx,3);
	  proC << ploty << "  ";
	}
      if (k==9)
	{
	  double plotx=fB1.answer[k].end-6;
	  double ploty=coe_0+coe_1*plotx+coe_2*plotx*plotx+coe_3*pow(plotx,3);;
	  proC << ploty << "  " << endl;
	}
      for(int s=0;s<100;s++)
	{
	  double plotx=fB1.answer[k].begin-6+s*0.01;
	  double ploty=coe_0+coe_1*plotx+coe_2*plotx*plotx+coe_3*pow(plotx,3);
	  double real=1/(1+plotx*plotx);
	  error1.push_back(abs(ploty-real));
	}
	
    }


  proC <<"----"<<endl;


  Bspline fB4(11,x1,y1);
  double f5_;
  f5_=(26.0*148)/pow(26,4);
  fB4.answer=fB4.S_23(2,f5_,f5_);
  proC << "-----" << "cubic Bspline D2" << "-----" << endl;
  
  for (int k=0;k<10;k++)
    {
      double coe_0=216*fB4.answer[k].coe[3]+36*fB4.answer[k].coe[2]+6*fB4.answer[k].coe[1]+fB4.answer[k].coe[0];
      double coe_1=108*fB4.answer[k].coe[3]+12*fB4.answer[k].coe[2]+fB4.answer[k].coe[1];
      double coe_2=18*fB4.answer[k].coe[3]+fB4.answer[k].coe[2];
      double coe_3=fB4.answer[k].coe[3];

      /*
      double x=fB4.answer[k].begin-6;
      double y=1/(1+x*x);
      cout << y <<" ";
      double yy=coe_0+coe_1*x+coe_2*x*x+coe_3*pow(x,3);
      cout << "** " << yy << " "<<endl;
      */

      for(int i=0;i<7;i++)
	{
	  if (fB4.answer[k].begin-6<xd[i] && xd[i]<=fB4.answer[k].end-6)
	    {
	      SD2[i]=coe_0+coe_1*xd[i]+coe_2*xd[i]*xd[i]+coe_3*pow(xd[i],3);
	    }
	}
  
      for(int s=0;s<10;s++)
	{
	  double plotx=fB4.answer[k].begin-6+s*sub;
	  double ploty=coe_0+coe_1*plotx+coe_2*plotx*plotx+coe_3*pow(plotx,3);
	  proC << ploty << "  ";
	}
      if (k==9)
	{
	  double plotx=fB4.answer[k].end-6;
	  double ploty=coe_0+coe_1*plotx+coe_2*plotx*plotx+coe_3*pow(plotx,3);;
	  proC << ploty << "  " << endl;
	}
      for(int s=0;s<100;s++)
	{
	  double plotx=fB4.answer[k].begin-6+s*0.01;
	  double ploty=coe_0+coe_1*plotx+coe_2*plotx*plotx+coe_3*pow(plotx,3);
	  double real=1/(1+plotx*plotx);
	  error2.push_back(abs(ploty-real));
	}
	
    }


  proC <<"----"<<endl;



  Bspline fB5(11,x1,y1);
  fB5.answer=fB5.S_23(3,0,0);
  proC << "-----" << "cubic Bspline natural" << "-----" << endl;
  
  for (int k=0;k<10;k++)
    {
      double coe_0=216*fB5.answer[k].coe[3]+36*fB5.answer[k].coe[2]+6*fB5.answer[k].coe[1]+fB5.answer[k].coe[0];
      double coe_1=108*fB5.answer[k].coe[3]+12*fB5.answer[k].coe[2]+fB5.answer[k].coe[1];
      double coe_2=18*fB5.answer[k].coe[3]+fB5.answer[k].coe[2];
      double coe_3=fB5.answer[k].coe[3];

      /*
      double x=fB1.answer[k].begin-6;
      double y=1/(1+x*x);
      cout << y <<" ";
      double yy=coe_0+coe_1*x+coe_2*x*x+coe_3*pow(x,3);
      cout << "** " << yy << " "<<endl;
      */

      for(int i=0;i<7;i++)
	{
	  if (fB5.answer[k].begin-6<xd[i] && xd[i]<=fB5.answer[k].end-6)
	    {
	      SD3[i]=coe_0+coe_1*xd[i]+coe_2*xd[i]*xd[i]+coe_3*pow(xd[i],3);
	    }
	}
  
      for(int s=0;s<10;s++)
	{
	  double plotx=fB5.answer[k].begin-6+s*sub;
	  double ploty=coe_0+coe_1*plotx+coe_2*plotx*plotx+coe_3*pow(plotx,3);
	  proC << ploty << "  ";
	}
      if (k==9)
	{
	  double plotx=fB5.answer[k].end-6;
	  double ploty=coe_0+coe_1*plotx+coe_2*plotx*plotx+coe_3*pow(plotx,3);;
	  proC << ploty << "  " << endl;
	}
      for(int s=0;s<100;s++)
	{
	  double plotx=fB5.answer[k].begin-6+s*0.01;
	  double ploty=coe_0+coe_1*plotx+coe_2*plotx*plotx+coe_3*pow(plotx,3);
	  double real=1/(1+plotx*plotx);
	  error3.push_back(abs(ploty-real));
	}
	
    }


  proC <<"----"<<endl;

  


  vector<double> x2;
  vector<double> y2;
  double n2=10;
  for (int i=1;i<=10;i++)
    {
      double x=i-5.5;
      double y=1/(1+x*x);
      x2.push_back(i+0.5);
      y2.push_back(y);
    }

  Bspline fB2(10,x2,y2);
  vector<double> plot2;
  fB2.answer=fB2.S_12(1.0/26,1.0/26);
  proC << "-----" << "quadratic Bspline" << "-----" << endl;

  for (int k=0;k<=9;k++)
    {
      double coe_0=36*fB2.answer[k].coe[2]+6*fB2.answer[k].coe[1]+fB2.answer[k].coe[0];
      double coe_1=12*fB2.answer[k].coe[2]+fB2.answer[k].coe[1];
      double coe_2=fB2.answer[k].coe[2];

      /*

      double x=fB2.answer[k].begin-6;
      double y=1/(1+x*x);
      cout<< x <<" ";
      cout << y <<" ";
      double yy=coe_0+coe_1*x+coe_2*x*x;
      cout << "** " << yy << " "<<endl;
      */

      for(int i=0;i<7;i++)
	{
	  if (fB2.answer[k].begin-0.5-6<xd[i] && xd[i]<=fB2.answer[k].end-0.5-6)
	    {
	      SD4[i]=coe_0+coe_1*xd[i]+coe_2*xd[i]*xd[i];
	    }
	}


	  for(int s=0;s<10;s++)
	    {
	      double plotx=fB2.answer[k].begin-0.5-6+s*sub;
	      double ploty=coe_0+coe_1*plotx+coe_2*plotx*plotx;
	      proC << ploty << "  ";
	    }
	  
	  if(k==9)
	    {
	      double plotx=fB2.answer[k].end-0.5-6;
	      double ploty=coe_0+coe_1*plotx+coe_2*plotx*plotx;
	      proC << ploty << "  " << endl;
	    }
	  for(int s=0;s<100;s++)
	{
	  double plotx=fB2.answer[k].begin-0.5-6+s*0.01;
	  double ploty=coe_0+coe_1*plotx+coe_2*plotx*plotx;
	  double real=1/(1+plotx*plotx);
	  error4.push_back(abs(ploty-real));
	}
      
      
    }
      
  
  proC <<"----"<<endl;

  vector<double> x3;
  vector<double> y3;
  for (int i=-5;i<=5;i++)
    {
      x3.push_back(i);
      y3.push_back(1.0/(1+i*i));
    }

  Bspline fB3(11,x3,y3);
  vector<double> plot3;
  fB3.answer=fB3.S_01();
  proC << "-----" << "linear Bspline" << "-----" << endl;

  for (int k=0;k<10;k++)
    {
      

      /*

      double x=fB2.answer[k].begin-6;
      double y=1/(1+x*x);
      cout<< x <<" ";
      cout << y <<" ";
      double yy=coe_0+coe_1*x+coe_2*x*x;
      cout << "** " << yy << " "<<endl;
      */

      for(int i=0;i<7;i++)
	{
	  if (fB3.answer[k].begin<xd[i] && xd[i]<=fB3.answer[k].end)
	    {
	      SD5[i]=fB3.answer[k].coe[0]+fB3.answer[k].coe[1]*xd[i];
	    }
	}

      
      for(int s=0;s<10;s++)
	{
	  
	  double plotx=fB3.answer[k].begin+s*sub;
	  double ploty=fB3.answer[k].coe[0]+fB3.answer[k].coe[1]*plotx;
	  proC << ploty << "  ";
	}
      if (k==9)
	{
	  double plotx=fB3.answer[k].end;
	  double ploty=fB3.answer[k].coe[0]+fB3.answer[k].coe[1]*plotx;
	  proC << ploty << "  " << endl;
	}
      for(int s=0;s<100;s++)
	{
	  double plotx=fB3.answer[k].begin+s*0.01;
	  double ploty=fB3.answer[k].coe[0]+fB3.answer[k].coe[1]*plotx;
	  double real=1/(1+plotx*plotx);
	  error5.push_back(abs(ploty-real));
	}
      
    }
  proC <<"----"<<endl;

  proC.close();

  ofstream proD;
  proD.open("proD.txt", std::ios::out);

  for (int i=0;i<7;i++)
    {
      double e1=abs(SD1[i]-yd[i]);
      Es1.push_back(e1);
      double e2=abs(SD2[i]-yd[i]);
      Es2.push_back(e2);
      double e3=abs(SD3[i]-yd[i]);
      Es3.push_back(e3);
      double e4=abs(SD4[i]-yd[i]);
      Es4.push_back(e4);
      double e5=abs(SD5[i]-yd[i]);
      Es5.push_back(e5);
    }
  proD << "-----" << "cubic Bspline D1" << "-----" << endl;
  for (int i=0;i<7;i++)
    {
      proD << Es1[i] << " ";
    }
  proD << endl;
  proD << "---"<<"plot"<<"---"<< endl;
  for(int i=0;i<=1000;i++)
    {
      proD << error1[i] << " ";
    }
  proD << endl;

  proD << "-----" << "cubic Bspline D2" << "-----" << endl;
  for (int i=0;i<7;i++)
    {
      proD << Es2[i] << " ";
    }
  proD << endl;
    proD << "---"<<"plot"<<"---"<< endl;
  for(int i=0;i<=1000;i++)
    {
      proD << error2[i] << " ";
    }
  proD << endl;

  proD << "-----" << "cubic Bspline natural" << "-----" << endl;
  for (int i=0;i<7;i++)
    {
      proD << Es3[i] << " ";
    }
  proD << endl;
    proD << "---"<<"plot"<<"---"<< endl;
  for(int i=0;i<=1000;i++)
    {
      proD << error3[i] << " ";
    }
  proD << endl;

  proD << "-----" << "quadratic Bspline" << "-----" << endl;
  for (int i=0;i<7;i++)
    {
      proD << Es4[i] << " ";
    }
  proD << endl;
    proD << "---"<<"plot"<<"---"<< endl;
  for(int i=0;i<=1000;i++)
    {
      proD << error4[i] << " ";
    }
  proD << endl;

  proD << "-----" << "linear Bspline" << "-----" << endl;
  for (int i=0;i<7;i++)
    {
      proD << Es5[i] << " ";
    }
  proD << endl;
    proD << "---"<<"plot"<<"---"<< endl;
  for(int i=0;i<=1000;i++)
    {
      proD << error5[i] << " ";
    }
  proD << endl;

  proD.close();
  
  
}
